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ABSTRACT 



The gravitational potential is a key function involved in many astrophysical problems. Its evaluation inside continuous 
media from Newton's law is known to be challenging because of the diverging kernel l/|r — r'|. This difficulty is generally 
treated with avoidance techniques (e.g. multipole expansions, softening length) themselves not without drawbacks. In 
this article, we present a new path that basically fixes the point-mass singularity problem in systems with, at least, two 
dimensions. It consists of recasting the gravitational potential ^ in an equivalent integro-differential form. 



■ip{r) 



fir 



■dq.g^'Hir) 



where (51,52) is a pair of independent spatial variables (linear and/or angular), / is a known function, and H is an 
auxiliary scalar function. In contrast with f/;, this "hyperpotential" H is the convolution of the mass density with a 
finite amplitude kernel k. We show that closed-form expressions for k can be directly deduced from the potential of 
homogeneous sheets. We then give a few formulae appropriate to the Cartesian, cylindrical and spherical coordinate 
systems, including axial symmetry. The method is essentially not limited, either on the geometry of the source or on 
the distribution, and its implementation is straightforward. Several tests based upon simple quadrature/differentiation 
schemes are presented (the homogeneous rectangular sheet, cuboid and disk, the Maclaurin disk and a truncated Lane- 
Emden solution). Compared with a direct summation, the extra computational cost is low and the gain is real: no 
truncated series, no free parameter, and a relative accuracy better than 1% for typically 16 nodes per spatial direction 
using the most basic numerical schemes. 
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1. Introduction 

Gravitation plays an essential role in the evolution of most 
astrophysical systems, from aggregates and dusty planetary 
rings to rotating stars, supermassive black holes in active 
nuclei and galactic clusters ([Ha chisu 1986; Koz hanovll2004l : 
iColwell et al1[2006t lKin3[20ia: .Comito et al...20lll) . In the 
investigation of various dynamical problems and equilib- 
rium configurations from first integrals and energy equa- 
tions, the potential appears as a fundamental scalar func- 
tion. In continuous media, it is naturally accessible through 
an integral, namely 



is difficult to estimate precisely the infiuence of any small 
element upon itself (i.e., "self-gravity"), which is not ame- 
liorate a lot by increasing the resolution. 

The point mass singularity can be avoided in vari- 
ous ways that are more or less faithful to Newton's law. 
The multipole expansion of the kernel \r — r'\~^ one of 
the most valuable theoretical tools in p otenti al theory 
CK ellogd [T929I: iDurandl [19531 iCohl et al.i [200lt lAksenovl 



ip{r) = -Q 



dm' 



\r — r' 



(1) 



where the kernel sweeps aways the point mass singularity 
— a direct consequence of Newton's inverse square law. 
Indeed, the potential integral is con vergent for niost den- 
sity distributions of physical inte rest ()Kellogaill929l : IDurandl 
Il953t iBinnev fc T remainel Il987l) . The singularity problem 
is inherent in the discretization-counting technique usually 
adopted. By dividing the system into small massive ele- 
ments and summing over all individual contributions, it 
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il999[ ). It is extremely efficient outside the material do- 
main and a few terms often suffice to reach computer pre- 
cision. Inside and even in close neighborhood, however, 
the convergence of the series is known to be poor be- 
cause r/r' w 1. Because the series is an alternate series, 
convergence is much delayed and truncations are critical 
(|Clementlll974l ). Low convergence is a common property of 
multipole expansions and is observed in various contexts 
other than gravitation (Wu ensche 1975; Kosov & Popelieg 
I2OOOI: iGramada fc Bournd I2OIOO . Users of muhipole ex- 
pansions generally need to incorporate a large number of 
terms — tens to h undred typically — before accuracy be- 
comes ac ceptable ("H achisul Il986t IStone fc NormanI 119921 : 
iMach fc Malec 2012). Because the number of integrals to 
estimate is equal to the number of terms, the computa- 
tional time increases linearly. Another option to derive 
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■0 is the Poisson 6' 

specifi c algo rithms ("Stone 
Spot^ 119951: [Briggs et al 



uation, which is ra 



Eidl:jL 



solved with 



Norman"1992' 'Storzer"199 
2000; Matsunioto & Hanaw; 



200.11: 1.TuseliusfcSundholmll2007HGuillet fc TevssieJl20l¥ 



Nevertheless, Poisson-solvers are not always "self-starting" , 
meaning they require precise boundary conditions only the 
integral approach can furnish. Another drawback is the 
shape of astrophysical bodies, which are often complex and 
not systematically match the n umerical meshes (for tech- 
niques based on mapping, see iGrandclement et al1l200lt 
iReesd 12006) . In particular, the Poisson equation is three- 
dimensional by nature, and not well-suited to problems in 
one and two dimensions. 

In this paper, we describe a novel path for determin- 
ing the Newtonian potential of a continuous system by re- 
casting Eq.([T|). The new form does not involve any singu- 
lar kernel, series, or "softening length" , but just the cross- 
derivative of the mass density convolved with a finite ampli- 
tude kernel. The recasting is exact and general in the sense 
that i) it preserves the Newtonian character of the interac- 
tions at all scales, and ii) it applies to any density distri- 
bution and morphology (shape and number of dimensions 
larger than one). This is therefore a new tool for both nu- 
merical applications and theoretical investigations in vari- 
ous domains of astrophyscis (e.g., simulations, generation of 
approximations, determination of potential/density pairs) 
and Physics as well. This paper goes bey ond the analy- 
sis presented in iHure &: Dieckmannl (J2012f) , which was re- 
stricted to axial symmetry. We consider here i) a generic 
treatment of the regularization step, regardless of the sys- 
tem of coordinates, ii) a full three-dimensional approach, 
iii) a simple recipe to determine the finite amplitude ker- 
nel, and iv) a direct application to the Cartesian, cylindri- 
cal, and spherical coordinates. 

The paper is organized as follows. In Section [5J we re- 
call the integral expression for the Newtonian potential of 
a continuous distribution. We formally describe the recast- 
ing of the potential integral based upon the properties of 
Newton's law (symmetry and independant spaces). The ap- 
plication to the Cartesian, cylindrical, and Spherical coordi- 
nate systems is the aim of Section [3] In particular, we illus- 
trate the method by considering a few test-cases, mostly 
of astrophysical interest, by using deliberately low-order 
numerical schemes (our goal is not to perform a critical 
study of the most efficient techniques for quadratures and 
differentiations). A conclusion summarizes the results and 
mentions possible issues to consider next. A few appendices 
contain formulae and demonstrations. 



P'Cq'i.Qj.qp 




boundary 



Fig. 1. Field point P(r), source point P'(r') belonging to 
the material domain J7', and elementary mass dm' . Inside 
ri', the separation PP' vanish. 



2.1. Idea behind recasting and strategy 

As Eq.(IT]) shows, there are two spaces in potential theory: 
i) the space of field points where the potential is requested 
(hereafter, the P-space), and ii) the space of source points 
that describes the source (hereafter, the P'-space). These 
spaces are superimposed in practice — this is the physical 
space — , but are decoupled mathematically. Indeed, when 
estimating the potential from Eq.((T]), r is held fixed while 
the integration is performed in the P'-space. The point- 
mass singularity, of hyperbolic-type, can be regularized us- 
ing two successive integrations in orthogonal directions (as a 
proof, note that the potential of flat or curved homogeneous 
sheets is generally a finite function of space and source pa- 
rameters). The idea is then to integrate the Newton kernel 
in the P-space until the singularity is finally suppressed. 
This operation is necessarily possible since the Newton ker- 
nel is symmetrical with respect to r and r' . Concretely, if 
P has coordinates {qi,q2,<l3), then the regularization step 
becometQ 



f{r)dqidq2 

PP'{qi,q2,q3) 



^9192 



(r;r') 



(2) 



where / is introduced for convenience (see below) and it 
is a function of P only. At this stage, the coordinates 
{qi',q2',q3') of P' are regarded as parameters, and k'^''^ 
must be a function of r. Since the regularization is per- 
formed in the P-space, it is made regardless of the mass 
distribution, which is especially attractive. The new kernel 
^9192 (hereafter the "hyperkernel" ) has, by construction, a 
finite amplitude and can be convolved with the mass den- 
sity. This is the convolution step: 



2. Recasting of the potential integral 

The Newtonian potential at a point P(r) in space of a body 
is given by Eq.([T|). The integral extends over the material 
domain il' (including its boundary), i.e., ip{r) = ■il;{r;il'), 
dm'{r') is the elementary mass at P'(r') e 51', |r— r'| =PP', 
and Q is the constant of gravity. The configuration is illus- 
trated by Fig.[T] As is well known, PP' vanishes everywhere 
inside J7', making the kernel singular, while the pot ential is, 
most of the time, a finite function of space (see e.g. iKelloggJ 
1192911 . 



K«i«^(r;r')dTO'EE-H(r;f7'), 



(3) 



where the factor —Q is introduced for convenience (see be- 
low) . This integral produces an auxiliary scalar function, % 
(hereafter, the "hyperpotential" ) . The Newtonian potential 
is then recovered by reversing the regularization-step. This 



^ Under invariance, the hyperbolic singularity can be con- 
verted into a logarithmic singularity that is subsequently regu- 



larized by a single integration, 
case (see Sect.[3|. 



I.e.. 



\r — r' I 



but this is a special 
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is the recovering step: 



dU^ = -5dU 



In 



n'i^'i'-dm' 



(4) 



^9192 



9192 



) dm' , 



dm' 



\r — r' 



= -Gf 

= f^- 

The advantage of this approach is twofold: the singularity 
is circumvented, and at the same time, it is accounted for 
exactly. In practice, the absence of diverging kernel renders 
step 2 easier than with the Newton kernel. Because convo- 
lutions produce smooth functions, step 3 is also expected 
to be uncomplicated. Step 1 is by far the most critical, but 
it is made once only provided the hyperkernel is analytical 
(there is no interest in the recasting if n'^^'^^ is to be deter- 
mined by numerical means). The gravitational potential is 
finally found from steps 2 and 3. The extra-cost is therefore 
low: there is only an additional differentiation compared to 
the classical approach, but the singularity is correctly man- 
aged. 
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coordinates of field points 
(P-space) 

field point P 



coordinates of source points 
(P'-space) 

source point P' 



{x,y,z) Cartesian coordinates {x',y',z') 

{R,9,z) cylindrical coordinates {R',d',z') 

{r,6,(f>) spherical coordinates [r' ,9' ,(f>') 



Fig. 2. Notations for the Cartesian, cylindrical, and spher- 
ical coordinates: P-space [left) and P'-space {right). 



2.2. Note on Chandrasekhar and Lebovitz superpotentials 

Our approa ch may evoke some aspects of the theory de- 
velopped in Chandrasekhar & Lebovitz ("1962) and subse- 
quent papers (see also Chandrasekhar 1987). These authors 
have shown that the Newtonian potential is the trace of a 
symmetric tensor potential, but is also the source term of 
a Poisson equation, namely 



V^x 



'2V, 



where the solution is 



X- 



I dm'. 



(5) 



(6) 



It is clear that H and x (called "superpotential" ) share two 
common properties: i) they are the convolution of the den- 
sity field by a finite amplitude kernel, and ii) they exactly 
reproduce the gravitational potential by partial differenc- 
ing. However, the present recasting differs from the theory 
of superpotentials on the following points: 

i) V is determined here through a single second-order par- 
tial derivative (not three); 

ii) the recasting is not limited to 3D-problems, but works 
for 2D-problems as well; 

iii) it is not really specific to a particular system of coor- 
dinates, while Cartesian coordinat es are mostly used in 
[Chandrasekhar fc Lebovita (|l962f ). 

In some sense, the hyperpotential H is some kind of op- 
timized version of the x-function, especially designed for 
numerical applications which is, initially, our main motiva- 
tion. 



3. Results and examples 

3.1. Derivation of hyperkernels, and the link with the 
potential of homogeneous sheets 

According to Eq. ([2]) , the recasting depends on the capabil- 
ity to determine analytically an expression for k^^"^^ asso- 
ciated with a given pair (91,^2) of coordinates, preferably 



a closed-form. There are many possibilities, in particular 
because of the presence of the function /(qi, 52, 93)7 which 
adds a degree of freedom. If we define / such that 



d'^A^ f{qi,q2,q3)dqidq2 



is an area element then 



^9192 



d^AjP) 
s \r-r'\'- 



(7) 



(8) 



where 5 is a surface 93 = const in the P-space. We see that 
Eq.® is nothing but, up to a factor —Q, the formula for 
the gravitational potential of a homogeneous surface with 
unit surface density dm/d^A, except that the role of the 
P-space and P'-space is exchanged. To obtain a formula for 
^9192^ it is sufficient to extract from the list of known poten- 
tial/density pairs those that correspond to a bi-dimensional 
distribution (i.e., a sheet) and constant surface density. As 
Eq.® suggests, there is no special constraint on the shape 
and size of the sheet (i.e., flat, curved, rectangular, circu- 
lar, etc.). Nevertheless, from a practical point of view, it 
seems preferable that S and fi' be geometrically "compati- 
ble" enough to facilitate the convolution-step. Since there is 
certain freedom in selecting the integral bounds in Eq.®, 
it is also better to consider a finite size (and finite mass) 
sheet. 



3.2. An easy implementation 

For each point P where ip is requested, the sequence of 
operations is the following: 

1. computing the hyperkernel k''^''^ from an appropriate 
formula (this depends on the coordinate system; see 
next section); 

2. estimating the hyperpotential H from Eqs.Q for the 
actual density distribution (a volume density p or a sur- 
face density E). A quadrature scheme is needed. This is 
the first place where numerical errors are generated; 
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qi,q2 f d^A surfaced hyperkernel K^^'^ 



Cartesian coordinates 


x,y 


1 


dxdy 


rectangular sheet 


K^y. see Eq.(|A.4|| 




y,z 


1 


dydz 


rectangular sheet 


K^" (see k"'') 




x,z 


1 


dxdz 


rectangular sheet 


k"" (see ti^y) 


Cylindrical coordinates 


e,z 


R 


Rdddz 


piece of hollow cylinder 


— 




R,e 


R 


RdOdR 


polar sector 


t.'^^seeY^q.]^ 




R,z 


1 


dRdz 


meridional sheet 


K^^, see Eq.(|B.2|) 


Spherical coordinates 


0A 


r^ sin <f) 


r^ sin <f)d<f)d6 


spherical cap 


— 




r,4> 


r 


rdrd(f) 


meridional sector 


k''*', see the Appendix |Dj 




r,e 


r sin (j> 


r sin 4>drd9 


piece of cone 




Axial symmetry* 


R 


R 


dR 


disk 


K^. see Ea.(|C.6|) 




r 


1 


dr 


cone 


— 







1 


d(j> 


spherical cap 


" 



Table 1. Pairs {qi, q-z), function /, and associated area element (PA for the three most popular coordinate systems. The 
formula for the hyperkernel k'^^'^^, when known in closed form, is indicated in the last column (otherwise ' — '). *Under 
axial symmetry, a single variable is necessary (see Sect. 




Fig. 3. The hyperkernel k^^ in the {x, y)-plane in the vicin- 
ity of point P' (magnified by a factor 5), for z = z' and 
three different values oi y — y' labelled on the curves. The 
Newton kernel which diverges for |r — r'| = is shown for 
comparison. 




-2.5 



-3.5 



Fig. 4. Error index e on potential values in the plane of 
the homogeneous square sheet with vertices at (±i, ±i, 0) 
(boundary in white). The mean error index is —2.52 
(dashed line). 



estimating the cross-derivative oi H. A differenti- 
ation scheme is needed. This requires determining 
hyperpotential- values in the vicinity of the actual point 
P. This is the second place where numerical errors are 
generated. 



Clearly, all techniques curently used to compute ip from 
Eq.([T]) can be employed for H. Because H is a convolu- 
tion (see below), fast specific algorithms coupled with high- 
perfomance differentiation schemes can probably be envis- 
aged. Tests presented in the following will employ the most 
basic schemes, which have produced good results. 



3.3. Results. The case of Cartesian coordinates 

Table[T]lists pairs (<Zi, <Z2), the associated area element (PA, 
and the function / appropriate for the Cartesian, cylin- 
drical and spherical coordinates that are most often used. 
Other coordinate systems and geometries can obviously be 
considered as well. The notations are summarized in Fig. 
[5] In Cartesian coordinates, the recasting that corresponds 
to the pair (51,52) = (2^,2/) is 



^(r) = dl. I K^ydm\ 



(9) 



and other pairs can be considered by permutation. As ar- 
gued above, we can determine k^^ by considering a surface 
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(73 = 2: = const, i.e., a flat horizontal sheet, and the most 
natural choice is the rectangular shape. The formula for 
the potential is known in that case (e.g. Durand 195^1). 
It is reproduced in Appendix |21 The hyperkernel is de- 
duced by exchanging P and P', which leads to Ea. (IA.4l) . 
Figure [3] displays k^^ in the {x, t/)-plane in the vicinity of 
the point P'{x' ,y'), for three different values oi y — y' and 
for z — z' = 0. The Newton kernel l/|r — r'| is also shown for 
comparison. We can see that the hyperkernel is a smooth 
function with finite amplitude. In particular, it is zero at 
zero relative separation. Any type of mass density profile 
can be injected in Eq.®, bi- or three-dimensional, uniform 
or not. 

We illustrate the potential recasting with two sim- 
ple examples. As a first test, we consider a square sheet 
in the {x, j/)-plane with constant surface density Sq (or 
dm' = T,odx'dy'), length unity, and centered on the origin, 
with vertices at (±i, ±i, 0). It is discretized on a A''' x M' 
grid with regular spacing in each direction. The P-grid 
is made of A^ x M points with uniform spacing as well. 
The quadratures and partial derivatives are determined by 
seco nd-order schemes that are among the most basic ones 
(e.g. iPress et al.lll992f) . Figure U shows the error index 
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Fig. 5. Same conditions and same color code as for 
Fig. m but for the homogeneous cuboid with vertices at 
(±i,±i,±i). The mean value is —3.00 (dashed line). 



log 



2 X 10^ 



-16 



1- — 



(10) 



between ip determined from Eq.(j9|) and the reference po- 
tential ipe (2 X 10~^^ is for double-precision computations). 
This case is for N' ^ M' = 32 and A^ = Af = A^' -^- 2, which 
leaves just one point outside the sheet, left and right, bot- 
tom and top. We see that the relative error is rather uni- 
form inside the material domain, of about 0.1%. Close to 
the edges of the sheet, the error rises slightly. Outside fl', 
the accuracy is still uniform, but typically better by an or- 
der of magnitude. 

As a second test, we consider a cube with uniform den- 
sity po (i-e., dm' = podx'dy'dz), length unity, and vertices 
at (±i, ±i, ±i). The numerical setup is the same as for the 
sheet, and z = 0. The reference potential is also known for 
this 3D body (|MacMillanlll"930t lWaldvogelllT976l) . FigureE] 
displays the error index versus x and y in the cube's mid- 
plane. Again, we notice that the relative error is uniform, 
with 0.2% typically inside the body, and a factor 10 better 
outside. Edge effects are less marked than in 2D. The in- 
tegration of the kernel in the third direction smoothes the 
errors, and the potential is now derivable when crossing the 
lateral faces of the cube. The Fortran 90 program used in 
these two examples is available upon request. 



3.4. The hyperpotential is a convolution product 

As shown in AppendixIXl when we set X = x—x', Y 
and Z = z — z' , the hyperkernel k^^ becomes 



y-y 



— Zatan 



XY 



Z\r — r' 



rin 



X 



a: In 



Y 



\r — r 



VX^ + Z2 



= K-yix,Y,z), 



where |r — r'| = vA^ + Y^ + Z^, and so the hyperpotential 
writes in the 3D case 



H{x,y,z) 



p{x\y',z') 



X K^^(x — x',y — y', z — z')dx'dy'dz' . (11) 



Since p = outside ft' , the integral bounds can be safely 
changed for ±00. We then conclude that H is a convo- 
lution product. This is expecte d because the potential it- 
self is a convolution product (jBinnev fc Tremaing 119871 : 
iHackbusch et al.ll2010[ ). We have 



d^,y{p^^^y) = p^dlyK^y 

1 

= p * 



(12) 



\r — r'\ 

This result is independent of the coordinate system, namely: 
E*k9i9^ in 2D, 



•H = 



(13) 



^p*K«i9% in 3D. 



3.5. Results in Cylindrical and Spherical coordinates 

In curved geometries, there are apparently fewer options. 
The reason is that a few formula for the potential of canon- 
ical surfaces are missing yet, and it is hard to find closed- 
form expressions for the hyperkernel by direct integration 
of Eq. ^ . One can probably use a series representation in- 
stead, but any truncation is expected to produce an approx- 
imate potential. Surfaces (73 = const of particular interest 
are (see Fig. [2] and Tab. [l]) 

1. in cylindrical coordinates: 

(a) a piece of a hollow cylinder (surface R = const); 

(b) a meridional sheet (surface 6 — const); 

(c) a polar sector (surface z = 0); 
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midplane (z-z'=0) z=z'=0.05 

1 / 





Maclaurin disk 
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0.2 0.4 0.6 0.8 1 



0.5 



Fig. 7. Same legend and same color code as for Fig. [SJ but 
a radially inhomogeneous sphere (boundary in white; see 
text). The mean index is —2.83 (dashed line). The P-grid 
{R, z) is a square, larger than the sphere's radius by 10%. 



Fig. 6. Error index for the homogeneous disk (top panel), 
and for the (inhomogeneous) Maclaurin disk (bottom panel). 
The set-up is the same in both cases (edge at R' — 1), and 
N' ~ 32. The computational grid has A^ = 33 points, but 
the same spacing. 



2. in spherical coordinates: 

(a) a piece of spherical shell (surface r — const); 

(b) a meridional sector (surface 6 = const); 

(c) a piece of cone (surface (f> = const). 

For cases la, 2a, and 2c (with </>' < |), the potential is 
apparently not known in closed-form; this would be helpful. 
We have no hyperkernel to propose. This question remains 
open. Case lb is accessible since the meridional sheet is 
nothing but a rectangular sheet in the plane {x, z), rotated 
counter-clockwise by an angle 9'. The formula for k^'^ can 
then be deduced from the Cartesian case (see Appendix IB]) : 
this is Ea. (jB.2p . Case Ic can also be treat ed since the p oten- 
tial of a polar sector has been derived in iHurel (|2012f ) . The 
formula is reproduced in the Appendix [Cl and the hyperk- 
ernel K^^ is given by Ea. (jC.4l) . Finally, case 2b is feasible 
since the meridional sector is a polar sector. We can then 
use the result for k^^ established in cylindrical coordinates 
and apply convenient rotations to derive k'"'^. A more direct 
calculus is presented in Appendix |D] 

3. 6. Axial symmetry 

If the source is axially symmetrical (i.e., de'P = 0), the 
integration of the Newton kernel in the P'-space over the 
polar angle 6*' e [0, 2tt] leads to 



de' 



\r — r' 



.K(fc), 



(14) 



where K is the complete elliptic integral of the first kind, 
^2 = (i?' + Rf + {z- z')^ and M = 2^/RR' . Because this 



function is hyperbolically singular when fc — )■ 1, we can de- 
termine an hyperpotential by an integration in the P-space. 
Again, there are not many options because we lack formulae 
for the potential of the hollow cylinder, for the cone, and for 
the piece of spherical shell (see Tab.[T]). Fortunately, there 
is the formula for the potentia l of the circular disk, i.e., a 
closed-form for f ^K(k)R ' dR' (iDurandl fToSSl: iKroueh et al] 
ll982l:lLass fc Blitzeill983tlHurel2012[ ) We can therefore de- 
duce an axially symmetrical hyperkernel k^ by exchanging 
the role of P and P' (see also Appendix [Ujl . This leads to 
Ea. (IC.6|) . The potential is also axially symmetrical, and it 
finally writes 

^ = ^dR-H, (15) 



R 



where 



n 



K dm' . 



(16) 



We now present three last examples. Figure El is for the 
circular disk with radius unity. The P-grid and the P'-grid 
are made of N' = 32 and iV = A^' + 1 points equally spaced 
in i?2 anfi jj'2 respectively — this radial scale is natural 
in this type of problem, both for the convolution and for 
the derivative. As above, the two grids coincide inside il' 
(there is just one point outside the disk) and the numer- 
ical schemes are second-order. The error index is shown 
at four different altitudes, including for the disk midplane 
(i.e. z = z'). The top panel shows the homogeneous disk. 
We see that the relative error is, on average, of about 0.1%, 
with again a slight degradation near the edge. The bot- 
tom panel shows the Maclaurin disk where E = \/\ — i?'^. 
For t his inhomogen eous case, the reference solution is taken 
from ISchuld (|2009[ ). We see that the relative error is about 
5 X 10"'^, which is not as good as in the homogeneous 
case. The error is almost insensitive to the altitude from 
the disk plane. This is due to the actual surface density 
profile (quadrature schemes are generally not very efficient 
to manage such situations). 



J.-M. Hure: Exact, singularity-free recasting of the Newtonian potential 



The last example is an inhomogeneous sphere with ra- 
dius unity. The potential/density pair is the following 



r' < 1 



p{r') - 



n(7rr') 
ttt' 
AG 



[i(7rr) 



(17) 



r' > 1 : p(r') = 

V'e(r) = -^ 



4G 1 

7r r ' 



which corresponds to the solution of the Lane-Emden equa- 
tion with polytropic exponent 7 = 2 (or index n = 1), 
truncated at the first zero. The hyperpotential in spherical 
coordinates writes 



H = 



sin0'd0' / p{r')n^r''^dr'. 
"'0 



(18) 



where k^ is the cylindrical hyper kernel k,^. This is a typical 
example where the surface S (a disk) and the domain Vt' 
(concentric shells) are somewhat disconnected. Here, the 
sphere is discretized into N' x M' points equally spaced 
in the (r', (/)')-plane. Potential values are determined from 
Eq. p3|) in the {B? , 0)-plane at A^x M points equally spaced; 
the computational box is larger than the sphere's radius by 
10%. Figure [7] shows the error index for N' = M' = 32 
and TV = M = 32. We see that the deviation is remarkably 
homogeneous inside and outside the sphere. The relative 
error is about 0.2% inside and outside the material domain. 
This is comparable to the case of the cube. 

Centrally symmetrical configurations can also be 
treated by using an hyperkernel, but there is nothing re- 
ally new here (see Appendix IE)) . 

4. Summary and concluding remarks 

We demonstrated that the Newtonian potential of con- 
tinuous bodies can be determined from the partial cross- 
differentiation of the mass density convolved with a finite 
amplitude kernel (a hyperkernel), regardless of any coordi- 
nate system. The recasting of Newton's integral is free of 
singularity and exact, and it applies to any type of two- 
and three-dimensional systems. Provided the hyperkernel 
is analytical, the extra-cost with respect to direct estima- 
tions is weak or negligible: it is only N operations vs. iV^ 
in a grid with N points. It is much lower if the method 
is only used to generate boundary conditions (and coupled 
with Poisson-solvers based on FFTs). The gain in accu- 
racy is huge since i) direct estimates cannot avoid errors, 
ii) there is no free parameter, and iii) there is no trunctaed 
series. We have given a few examples in Cartesian, cylindri- 
cal, and spherical coordinates that prove the efficiency even 
with low-order quadrature and differentiation schemes. The 
recasting is therefore very attractive for any numerical ap- 
plications. It is also a new tool for investigating various 
theoretical problems and derive potential/density pairs or 
approximations . 

This work can therefore be continued and improved in 
several ways. Knowing analytical expressions for the hy- 
perkernel associated with a given pair of orthogonal coor- 
dinates is the critical point of the method. We have shown 
that hyperkernels can be directly generated by consider- 
ing the potential of homogeneous sheets, while there are 
doutblessly other possibilities. For Cartesian coordinates, 



all hyperkernels are known. In cylindrical and spherical ge- 
ometries, a few closed- form expressions are apparently lack- 
ing yet (hollow cylinder and cone for instance). It would 
therefore be interesting to investigate this kind of question. 
The formula for the polar sector should be helpful for most 
astrophysical applications however, such as for modelling 
rotating fluids. Other coordinate systems and geometries 
can be envisaged. For instance, the potential of inhomoge- 
neous elliptic bodies can be determined, as done under cen- 
tral symmetry (see the end of Sect. 13. 6p , from the theory of 
thin homeoids (e.g. lMacMillan|[l930t ) . This show the impor- 
tance of seeking new potential/density pairs associated to 
2D-systems. It would also be interesting to analyze in more 
detail the numerical implementation of the method. There 
is obviously a wide panel of techniques at our disposal to 
perform quadratures/convolutions and differentiations, fi- 
nite differences (as considered here), spectral methods, etc. 
Finally, applications exceed the astrophysical context 
of gravitation. The approach is obviously suited to electro- 
statics and to electromagnetism since the potential vector 



Air) 



udq' 



\r — r' 



(19) 



where u is the velocity of electric charges. It is also trans- 
posable to incompressible hydrodynamics where the pres- 
sure p obeys a Poisson equation too. 



p{r) = - 



V • [(m • V) u] dm' 



(20) 



where u is the fluid velocity. 
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Appendix A: Hyperkernel for the rectangular sheet (Cartesian coordinates) 

The notations are those of Fig.[51 Up to a factor —Q, the potential of a homogeneous rectangular sheet with unity surface 
density is found from the integral 

dx'dy' ^^^^^ 



y/{x' - xy + [y' - yY + {z' - z)2 

where the bounds represents the coordinates of the four corners of the rectangular sheet. A closed- form expression is 
found for instance in lDurandl ()1953[ ). If we set X = x — x\Y = y — y' and Z = Z— z', we have |r — r'| = \/ X'^ + Y'^ + Z'^ 
and dx'dy' = dXdY. The indefinite integral is 



'^'^'^^^ =X-Y\n{\r-r'\-X)-X In (|r -r'\-Y)-Z 



\r — r' 



X XY 

atan — - + atan 



Z Zr-r' 



(A.2) 



For the present prupose, we only need to generate the Newton kernel by a mixed partial derivative d^ /dxdy, and so, we 
have a certain liberty in choosing the most convenient integral bounds. Here, we take x and x' for the integral over a;', 
and y and y' for the integral over y. from this, we obtain 

"^'"^ f^'^y dY \r-r'\-X \r-r'\-Y XY 

dX -^^ = -Vln l'^ "^1 ^ -yin '"^ "^1 _L _Z.t..n (A.3) 







\r-r'\ ^Y"^ + Z2 ^pG + Z^ Z\r - r' 



Consequently, the hyperkernel is obtained from this expression by exchanging the variables x' and x, and y' and y. We 
find 

\r-r'\+X \r-r'\+Y XY 

K^^^rin^ ' +Xln '\ ' ^-Zatan , , . (A.4) 

VF2TZ2 ^/WTZ^ Z\r-r'\ ^ ' 

To generate k^^ associated with a potential expressed as d'^T-L/dydz, we perform the permutations {x, x') <-^ {z, z') in 
the formulae above. To generate k^^ associated with a potential expressed as d'^'H/dxdz, we perform the permutations 
iy,y')^ iz,z'). 

Appendix B: Hyperkernel for the meridional sheet (cylindrical coordinates) 

The notations are those of Fig. [2j First, we use Ea. (|A.4p and perform the permutations {y,y') O {z,z'). We find 

n-^Zlnt^Z:i±£+Xlnt^i:i±£-Y.t.n-^. (B.l) 

VY^ + Z2 V^2 + y2 Y\r-r'\ ' 

Then, we apply a counter-clockwise rotation (i.e., positive trigonometric sense) around the 2:-axis by an angle 9'. We 
obtain the meridional sheet. The hyperkernel then writes 



= (z - z) In , - - {R' + R cos 2;9) In 



^R^sm'2(3+(z~z')^ VR' + R"-2'R'Rcos2l3 

r. ■ c.n {z - z')(R' + Rcos2P) ,^^, 

-l-i? sin 2/3 atan -i '--^—, — ^, B.2 

^ i?sn2/3 r-r' ' ^ ' 



and the potential is given by 

" ' "" ~ ■ '■'Rz 



^(r) = -Gdl, / K^'^dm'. (B.3) 



O' 



Appendix C: Hyperkernel for the polar sector (cylindrical coordinates) 

The notations are those of Fig. [51 Again up to a factor —Q, the potential of a circular sector is found from the formula 

R'dR'de' 



|r — r' 

where 

„/|2 _ / D' I D^2 I / /^2 



\r - rf = {R' + RY + [z - z'f - AR'Rcos^ 



2 
This double integral has been calculated in iHurel ()2012l) . The indefinite form is 

p/^p/^iQ/ p/2 p2 p/ /? /*2 

I ,1 = 5E{(3,k) + ^ F(/?,fc) + :g_g^n(/3,m2,fc) - i? sin 2/3 asinh 

C(i?' + i?cos2/3) 



^atan 



i?sin2^ \r-r'\ 





(C.l) 




(C.2) 


R' + Rcos2(3 




^/C + R^ sin^ 2/3 






(C.3) 
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where F{ip, k), E{(p, k) and !!((/), to^, k) are the incomplete elliptic integral of the first, second, and third kinds, respectively, 
5={R' + RY + C'^,C = z- z',k5 = 2\/R'R, 2/3 = n - {9 - 9'). To generate k™, it is sufficient to consider the following 
integral bounds: and R' for the radial integration, and 9 — tt and 9' for the angular part. Next, the source point and 
the field point are exchanged (note that S, k, ^^ and m are not impacted). We finally obtain 



K 



d2 d'2 d p/ /-2 

5 ^'^' '^ R + R' 6 



6E{f3', k) + F(/3', k) + -—— :yn(/3', m^, k) (C.4) 



r., ■ ^a • , R + R'cos2fi . , R'cos2fi 

+ R sm 2p asmh — ^^^^^^^^^= — asmh — r^^^^=^^^^= 
\ yJC + ^'2 sin' 2/3 v/C^ + i?'2sin'2/3^ 



C < atan 



C(i? + i?' cos 2/3) 



a sin 2/3 \r — r' 
where 2/3' — t: — {9' — 9) —2tt — 2/3. The potential is then given by 



/ c ^ 

atan — - pntan 2/3 



Hr)=-^gdlel n^'dm'. (C.5) 



Under axial symmetry, we have 



/.^ 



MCV + ^E(fc) + :^l-^K(fc) + ^§_|n(m^ fc) 



(C.6) 



(5 ^ ' 5 R + R' 

where E, K and 11 are the complete elliptic integrals of the first, second, and third kinds, respectively. 

Appendix D: Hyperkernel for the meridional sector (spherical coordinates) 

The notations are those of Fig. [21 Up to a factor —Q, the potential of a meridional sector defined by 9' = const, is found 
from the double integral 



// 



r — r' 



with convenient bounds. We can calculate this expression directly from the formula for the polar sector in cylindrical 
coordinates (see Appendix ICj) . First, we write the relative separation |r — r'| in the following form 

|r-r'p = [r' + rvf +r'^{l-v'^)-'ir'rvsin^T, (D.2) 

where 

V= Vl - sin' </>' sin' 2/3 

2t = TT - (00 - <P) 

tan 00 = ~ tan 0' cos 2/3. 
Then, we notice that the similarity between Eq. ljC.ip and Ea. (|D.l[) is perfect if we make the following substitutions 

■R' or' 
R^ rv 

{z- z'f Or'(l-i/') 

,/3 0T, 

and then 

f(5'=r'+r'' + 2rr'iy, 

i 1,2 _ 4rrV ^2 _ 4rr't/ (-'-'•'^j 



'^2 . . .2,1 „2^ (D.4) 



From Ea. (jC.3l) . we see that Eg. ljD.ip becomes 

r'dr'dS' ^ , ,, r'" — (rvY , ,^ r' — ri^ r'(l — j/')^ , , ,, , r' + rz^cos2r 

^ ^6E{T,k) + i-^-i^(T,fc) + !^ ^n(T,TO',fc)-ri^sn2rasnh 

|r-r'| ^ ' ^ 5 ^ ' ^ r' + ri. 5 I > ' ^ rVl - v'' cos' 2t 

^ rVT^^atan V T^C^ + r. cos 2r) ^ 
z^sm2T |r — r'| 

We then derive k''"^ ~ JJ I'j.Lj'i from Ea (jC.4p by making the same substitutions, and the potential is given by 



^{r)^--gdl^j K'-'f'dm'. (D.7) 

*" Jn' 
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Appendix E: Central symmetry 

When dr' p — 0, we can derive an hyperkernel by considering the potential of a spherical shell with unity surface density. 
From the Gauss theorem, we easily find (still up to —Q) 

sin 0'd0' / ^^ - ^H{r' - r) + —H{r - r'), (E.l) 

J2-K \i' - ■^ I r' r 

where H is the Heaviside function. The hyperkernel k^^^^^ is then obtained by exchanging r and r' in this expression. In 
this case, we have simply ip = H with 



■H = -g p{r'y^K''^''"dr' (E.2) 

J r' 

= -AttG j p{r')r'H{r' - r)dr' - — f p{r')r'^H{r - r')dr' , 

where the integral boun ds are the inner radius and outer radius of the sphere. This result is well known (e.g. 
iBinnev fc TremaindfToSTl) . 
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